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[i] A three-dimensional numerical simulation of plasma density irregularities in the 
postsunset equatorial F region ionosphere leading to equatorial spread F (ESF) is 
described. The simulation evolves under realistic background conditions including 
bottomside plasma shear flow and vertical current. It also incorporates C/NOFS satellite 
data which partially specify the forcing. A combination of generalized Rayleigh-Taylor 
instability (GRT) and collisional shear instability (CSI) produces growing waveforms with 
key features that agree with C/NOFS satellite and ALTAIR radar observations in the 
Pacific sector, including features such as gross morphology and rates of development. 

The transient response of CSI is consistent with the observation of bottomside waves with 
wavelengths close to 30 km, whereas the steady state behavior of the combined instability 
can account for the 1 00+ km wavelength waves that predominate in the F region. 
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1. Introduction 

[ 2 ] Around sunset, the loss of dominance of the E region 
dynamo to the F region creates a shear in the horizontal 
plasma flow at the F region bottomside [Haerendel et al., 
1992], Above the shear node, the vertical electric field 
points downward, and the plasma drifts eastward. Below 
the node, the vertical electric field points upward, and the 
plasma drifts westward, even though the neutral flow is 
everywhere eastward. The imperfectly efficient F region 
dynamo sets up a vertical current supplied by field-aligned 
currents connecting the equatorial F region to the E layer at 
low latitudes to midlatitudes. Related to this phenomenon, a 
family of plasma irregularities collectively referred to as 
equatorial spread F (ESF) disturbs the nighttime equatorial F 
region (see Woodman [2009] for details). Those irregulari- 
ties can be detected by incoherent and coherent scatter radars 
after sunset and appear in three different classes (or stages): 
topside, bottomside, and bottom-type. 

[ 3 ] Topside irregularities are driven by the generalized 
Rayleigh-Taylor instability (GRT) near and above the F 
peak. Typical GRT growth rates at the L 7 peak are ~ 10" 3 s -1 
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[Sultan, 1996; Hysell et al., 2006], which implies no more 
than 3-4 e-folds per hour. This is a universal instability 
(i.e., nearly scale independent), and thus GRT theory cannot 
explain by itself the predominant scale sizes of ESF 
irregularities. 

[ 4 ] At altitudes below the plasma shear node, the differ- 
ential motion between plasma and neutrals and associated 
vertical current provide a source of free energy that can drive 
collisional shear instability (CSI). During the initial transient 
phase, the zonal wavelength of the fastest growing mode (A) 
is approximately twice the density gradient scale length 
(La 1 = 1 1[), i.e., A — 2 L. Longer wavelengths dominate 
during the asymptotic phase, with A — 4ttL [Hysell and 
Kudeki, 2004]. For typical bottomside gradient scale 
lengths, the zonal wavelengths for the transient and asymp- 
totic growth phase are expected to be ~30 km and 100+ km, 
respectively. The irregularities produced by CSI propagate 
obliquely at an intermediary angle between upward and 
westward. CSI growth rates are faster than that of GRT 
under nominal conditions. Typical growth rates are ~5 x 
10 -3 s _1 , or 18 e-folds per hour during the transient phase 
[Kudeki et al., 2007]. CSI therefore can account for the scales 
and rapid growth of the dominant irregularities observed by 
radar in ESF. Bottom-type layers are almost ubiquitous over 
Jicamarca immediately after sunset, and their patchy signature 
seem to signal the onset of the CSI [Hysell et al., 2004]. 
Bottomside irregularities generated by CSI remain confined 
to the base of the F region, however, and are therefore not 
space weather effective. In order to bridge the gap between 
irregularities at the bottomside and the topside, zonal electric 
fields are generally needed to propel irregularities toward the F 
peak, leading to topside spread F. 
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[5] In recent years, efforts to simulate ESF in three 
dimensions have been made [Huba et al, 2008; Retterer, 
2010; Keskinen, 2010]. These simulations are driven in 
such a way as to emphasize the role of GRT while mainly 
neglecting CSI. More recent work has shown that the two 
instabilities function together in nature: CSI appears neither 
to suppress GRT nor merely to seed it; rather, the two 
mechanisms appear to grow and develop faster when acting 
together than either when acting alone, exhibiting a greater 
mixing depth [Aveiro and Hysell, 2010]. 

[6] The numerical simulations developed by Aveiro and 
Hysell [2010] incorporated a number of nonessential sim- 
plifications. Their model included time-invariant neutral 
atmosphere and background forcing. Also, the simulations 
were based on a single ion fluid model which did not take 
recombination into account. Here, we present the results of a 
multi-ion simulation which includes recombination, Coulomb 
collisions, a time-varying neutral atmosphere, and variable 
forcing, inferred from data from the C/NOFS satellite. The 
results are compared to incoherent scatter radar observations 
made by the ALTAIR radar in the Pacific sector. 

[ 7 ] Key attributes of our simulation included (1) its sen- 
sitivity to the background circulation, (2) its sensitivity to 
event sequencing, and (3) its inclusion of parallel electric 
fields. The significance of the latter is not evaluated explic- 
itly here and is left for a forthcoming analysis. However, it is 
worth mentioning that parallel ambipolar electric fields 
essentially changes the nature of the perpendicular fields, 
affecting convection and circulation [see, e.g., Zalesak etal, 
1990; Oppenheim and Dimant, 2006], The main issue we 
address below is the degree to which it is possible to predict, 
quantitatively, the behavior of ESF given an appropriate 
specification of the initial conditions and the forcing. 

2. Radar Observations 

[8] We begin with a presentation of the radar observations 
that will ultimately be used to motivate and validate the 
numerical simulations. The ARPA Long-Range Tracking 
and Instrumentation Radar (ALTAIR) is a VHF/UHF radar 
located at Kwajalein Atoll (8.8°N, 167. 5°E) [Tsunoda et al., 
1979]. Figure 1 displays ALTAIR UHF (422 MHz) data for 
the night of 29 April 2009. ALTAIR observations made 
perpendicular to the geomagnetic field at F region heights 
represent a combination of coherent and incoherent scatter, 
whereas off-perpendicular scans reflect only incoherent 
scatter. Some coherent scatter can be informative in data sets 
like these because it highlights fine structure that might 
otherwise be difficult to perceive. 

[ 9 ] Figure 1 (top) shows an ALTAIR scan made pointing 
perpendicular to the geomagnetic field around 09:45 UT 
(LT ~ UT + 1 1 h in the radar longitudinal sector). By this 
time, irregularities had formed at the base of the F region as 
bottom-type layers. These were patchy and spaced by about 
35 km. They are common precursors of fully developed 
ESF, whereas bottom-type layers tend to be more continuous 
(beam filling) or absent on nights when topside ESF does 
not develop [. Hysell et al., 2004]. 

[ 10 ] ALTAIR scans made pointing off perpendicular to 
the geomagnetic field at 10:30 UT and 11:26 UT are shown 
in Figures 1 (middle) and 1 (bottom), respectively. By about 
10:30 UT, the irregularities were still confined to the base of 


the F region, but the predominant horizontal scale had 
increased to 100-200 km. Between about 10:30 and 11:26 
UT, the bottomside depletions underwent strong vertical 
development. Depleted channels penetrated to the topside 
while retaining their 100-200 km horizontal spacing. The 
depletions had characteristic tilts and showed evidence of 
bifurcation and secondary instabilities. Wedges of enhanced 
plasma also appeared to be driven downward in between 
ascending depletions. 

3. Satellite Observations 

[ 11 ] To drive simulations of the 29 April 2009 ESF event, 
in situ electric field measurements from the Ion Velocity 
Meter (1VM) and Vector Electric Field Instrument (VEFI) 
on board the Communications/Navigation Outage Forecast- 
ing System (C/NOFS) satellite were used. This satellite was 
launched in April 2008, to investigate and forecast ambient 
plasma density and irregularities in the equatorial F region 
ionosphere [deLa Beaujardiere, 2004]. IVM is composed of 
an ion drift meter and an ion-retarding potential analyzer 
that together provide measurements of the ion drift vector, 
the ion temperature, and the major ion composition. The 
VEFI instrument consists of an electric field meter, a flux- 
gate magnetometer, an optical lightning detector, and a fixed 
bias Langmuir probe. Among its products are measurements 
of DC, quasi-DC, and AC electric fields. 

[ 12 ] Figure 2 shows the zonal electric field (based on IVM 
ion drift and VEFI magnetometer measurements) and (IVM) 
ion density measurements with the respective orbital infor- 
mation tabulated below (geographic and magnetic coordi- 
nates, altitude, and universal and solar local times). The 
results correspond to C/NOFS passes during the postsunset 
period on 29 April 2009, centered at the Kwajalein longitu- 
dinal sector. For comparison, satellite passes at —09:03 UT 
(Figure 2a), -10:47 UT (Figure 2b) and -12:31 UT 
(Figure 2c) are shown. 

[ 13 ] During the first C/NOFS pass at about 09:03 UT, the 
zonal electric field was essentially zero, and no density or 
electric field perturbations were detected in the altitude 
range between 495 and 520 km. By the time of the next pass 
at about 10:47 UT, the zonal electric field had increased to 
0.5-0.75 mV m _1 , and small electric field irregularities were 
emerging between 435 and 450 km altitude. That the back- 
ground plasma density increased as the satellite altitude 
decreased suggests that it was flying in the topside, which is 
consistent with the evidence of Figure 1. 

[ 14 ] By the time the satellite passed through the Kwajalein 
sector again at about 12:30 UT, the electric field irregulari- 
ties had become large (up to about 4 mV m” 1 ). Modest 
plasma density irregularities were also present at altitudes 
near 400 km. In view of Figure 1, however, the satellite 
most likely flew just above the most intense plasma den- 
sity irregularities. 

4. Numerical Simulations 

[ 15 ] Below, we describe a numerical scheme for the sim- 
ulation of the plasma density irregularities in the equatorial 
spread F event. The 3-D simulation updates the plasma 
density and electrostatic potential in time assuming quasi- 
neutrality and momentum balance. Inertia is neglected in the 
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Figure 1. (left) ALTAIR radar scans for 29 April 2009, showing backscatter power, range corrected and 
scaled to electron density on a logarithmic scale: (top) perpendicular scan at 09:45 UT and off-perpendicular 
scans at (middle) 10:30 UT and (bottom) 11:26 UT. (right) Vertical electron density profiles. Note that 
LT ~ UT + 1 1 h. 
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Figure 2. Ion density (cm ') and zonal electric field (mV m ' ) 
at three different passes through the Kwajalein sector: 
(a) -09:03 UT, (b) -10:47 UT, and (c) -12:30 UT. The tab- 
ulated entries are universal time (UT), magnetic latitude 
(MLAT), solar local tune (SLT), geographic latitude (GLAT), 
longitude (GLON), and altitude (ALT) in km. 


present incarnation. The characteristics of the neutral atmo- 
sphere are imported from empirical models. The background 
electric field is specified and controls the forcing. Simula- 
tions are performed with realistic background conditions 
including bottomside plasma shear flow and the attendant 
vertical current. 

[i6] Rather than assuming equipotential magnetic field 
lines, the simulation solves the potential equation in three 
dimensions, allowing us to keep track of the distribution of 
parallel currents forced by the dynamos and their closure. 
This also accommodates the possible effects of drift waves 
and related instabilities that may play secondary roles in ESF 
[Huba and Ossakow, 1979; LaBelle et al, 1986; Hysell 
et al, 2002]. For instance, the full 3-D treatment of warm 
plasma shows that perturbations on the surface of iono- 
spheric plasma clouds propagate and twist into a barber pole 
configuration [Drake and Huba , 1987; Drake et al., 1988; 
Zalesak et al., 1990]. As a consequence, the phase of the 
corrugations varies strongly along B, the plasma cloud loses 
its flute-like characteristic (i.e., k ■ B V 0), and the resulting 
finite A'n has a dissipative influence on the instability. We do 
not know the extent to which similar phenomenon might 
affect the ESF development, but a study is underway. 
Finally, it is worth mentioning that magnetic field lines cease 
to be equipotentials near the bottom of the E region, decou- 
pling the plasma and decreasing the conductive loading on 
the dynamo, which may influence the overall circulation. 

4.1. Conservation Laws 

[n] Neglecting stress, heat flow, and inertia, the fluid 
continuity equation and momentum equations for electrons 
and ions are given as 

-£+V-(n a v a )=P a -L a , (1) 


0 = - Vp a + n a q a { E + v„ x B) 


+ m a n a 


g - v an {\ a 


at/3 

U) — ^ u a p(\ a — , 

0 


( 2 ) 


where p a = n a k B T a is the pressure. The variables n a , q a , m a , 
T a , P a , L a , and v Q are number density, electric charge, mass, 
temperature, production and loss rates, and velocity of the 
species a, respectively. The loss term includes charge 
exchange and dissociative recombination for a given spe- 
cies. Nighttime photoionization is neglected, but production 
due to charge exchange is considered. Flere, /3 represents the 
other ionized bodies that exchange momentum or energy 
with a. The terms v a p and v a „ are the collision frequencies 
with charged particles and neutrals, respectively. The terms 
E, B, U, g, and k B represent the electric field, magnetic field, 
neutral wind, gravity, and the Boltzmann constant, respec- 
tively. Equation (1) is not in conservation form, since it 
includes chemical processes. However, it can be converted 
to conservative form through variable transformation. Two 
simplifications are applied to equations (1) and (2): (1) the 
neglect of diagmagnetic drifts (V« a x B), since their asso- 
ciated flux divergence is very small when compared to the 
other terms, and (2) Coulomb collisions (the last term on the 
right-hand side of equation (2)) are treated implicitly and 
only included in the parallel direction. 
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4.2. Numerical Scheme 

[is] The simulation algorithm performs two computations. 
First, the self-consistent electric field is found using a 3-D 
potential solver based on quasi-neutrality, which can be 
expressed as 

V ■ [S ■ V$] = V • [E - (E„ + U x B) + D ■ V« + S - g] , (3) 

where E 0 is the background electric field, E is the conduc- 
tivity tensor, D includes the diffusivity tensor, 5 • g is gravity- 
driven current density, and $ is the electrostatic potential 
that must arise to preserve a solenoidal current density. 
Explicit definitions of the tensors are given by Shume et al. 
[2005]. Note that the resulting electric field is the combi- 
nation of the background and electrostatic fields (i.e., E = 
E 0 — V$ ). Equation (3) involves a partial differential 
equation for the potential in three dimensions. It is solved 
using the BiConjugate Gradient Stabilized (BiCGStab) 
method [e.g., van der Vorst, 1992] using the algorithms 
described by Saad [1990]. 

[ 19 ] Second, we solve a discretized version of equation (1) 
for each ion species using a monotone upwind scheme for 
conservation laws (MUSCL) (a pedagogical review is given 
by Trac and Pen [2003]) directly applicable to the ion con- 
tinuity problem (once the recombination term has been 
absorbed into a redefinition of n a so that equation (1) 
becomes a conservation equation, as necessary). The char- 
acteristic of the neutral atmosphere (densities, temperature, 
and wind velocity) are updated in time on the basis of inputs 
from climatological models. The background electric field is 
specified and partly controls the forcing. The ion composi- 
tion includes 0 + , NO + , and 0 2 . The model was constructed 
using tilted magnetic dipole coordinates (p , q , (/)), where the 
tilt is matched to the magnetic declination in the longitude of 
interest. In our terminology, p represents the Mcllwain 
parameter (L), q is the magnetic colatitude, and 0 is longi- 
tude [see, e.g., Hysell et al., 2004]. 

4.3. Empirical Model Drivers 

[ 20 ] To initialize the model runs, we derive plasma num- 
ber densities from the Parametrized Ionospheric Model 
(PIM), a parametrization of the output of several regional 
theoretical model outputs generated for different climato- 
logical conditions and tuned somewhat to agree with data 
from Jicamarca and elsewhere [ Daniell et al., 1995]. PIM 
reproduces electron density profiles from Kwajalein and 
Jicamarca fairly closely if scaled slightly to account for day- 
to-day variability [Aveiro and Hysell, 2010]. The back- 
ground number density profiles used for initial conditions 
are intended to reproduce measurements from ALTAIR, 
including the F region peak height and peak density, 
the bottomside density scale height, and the density in the 
valley region. 

[ 21 ] Plasma mobilities are calculated using neutral com- 
position and temperature estimates from the Mass Spec- 
trometer and Incoherent Scatter (NRL-MSISE00) model 
[Picone et al, 2002] and ionospheric composition estimates 
from the International Reference Ionosphere (IRI-2007) 
model [Bilitza and Reinisch, 2008]. We take the ionosphere 
and neutral atmosphere to be in thermodynamic equilibrium 
after sunset (T„ = T e = T,). Expressions for the ion-neutral 


and electron-neutral collision frequencies used to compute 
conductivities can be taken from Richmond [1972]. 
Expressions for the Pedersen, Hall, and direct mobilities and 
diffusivities themselves are given, e.g., by Kelley [2009]. 

[ 22 ] For our runs, the zonal neutral winds are obtained 
from the Horizontal Wind Model (HWM-07) [Drob et al., 
2008]. As with PIM, these may be tuned slightly to 
account for day-to-day variability. The model is sensitive to 
the wind speeds throughout the thermosphere, since the 
growth rate of the collisional shear instability depends on the 
counter streaming flow of the neutrals and plasma. Meridi- 
onal and vertical winds were neglected in the simulations 
presented here. 

4.4. Simulation Setup 

[ 23 ] Our simulation was cast on a rectangular grid 139 x 
133 x 189 points wide in (p, q, cji) space constructed using a 
tilted magnetic dipole coordinate system. A cut through the 
equatorial plane spans altitudes between 90 and 510 km and 
longitudes of ± 6°. The flux tubes covered by the parallel 
coordinate all reach to the lower E region. The runs shown in 
this section are centered on the dip equator near Kwajalein 
(5.5°N latitude, 166. 5°E longitude, 7.4° declination) with 
background conditions modeled for 29 April 2009. To seed 
the simulation run, we added independent Gaussian white 
noise to the initial number density with a 20% relative 
amplitude. The local time dependence of the background 
zonal electric field (E 0 ) was estimated using measurements 
from IVM and VEFI on board the C/NOFS satellite. E, b 
was taken to be constant in the F region but its amplitude 
is reduced below the valley to avoid the growth of elec- 
trojet gradient drift instabilities. Based on the C/NOFS 
measurements E $ was modeled as 0 mV m _1 before 
21:00 LT, and 0.75 mV m~' after 22:00 LT with a gradual 
transition using a hyperbolic tangent function in between. 
For example, E^ = 34 pV m -1 on the western boundary and 
0.75 mV m -1 on the eastern boundary of the simulation 
by 10:40 UT, since those locations represent 21:22 LT and 
22:10 LT, respectively. 

5. Simulation Results 

[ 24 ] The initial conditions of our simulation at to = 
9:30 UT (=20:30 LT for the Kwajalein sector) are shown 
in Figure 3. Figure 3 (top) shows the background plasma 
density. The initial conditions are consistent with ALTAIR 
radar observations for 29 April 2009, at that time. Figure 3 
(middle) shows the perpendicular current density with the 
magnitude and direction indicated by the color disc in 
Figure 3 (bottom). There are two major components to the 
current density that are combined in this depiction. First, a 
zonally eastward current driven by gravity exists throughout 
the region, showing a maximum at the F peak near 330 km 
altitude. In addition, a zonally eastward A region wind drives 
a vertical upward current in the region above about 200 km. 
The empirically modeled wind is almost constant with alti- 
tude, and thus the vertical current scales like the background 
ion density, showing a peak near 330 km and a strong ver- 
tical gradient in the bottomside between 280 km and 240 km 
altitude. Below 150 km, the current is influenced by both the 
local zonal winds and the enhancement of the zonal current 
in the electrojet. 
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Figure 3. Initial conditions for ESF model run (i.e., at t 0 = 09:30 UT). (top) Plasma number density, 
(middle) Transverse current density in the equatorial plane with (left) equipotential contours (in kV) 
superimposed in white and (right) vertical electric field, (bottom) Current density in the meridional 
plane. The current densities are vector quantities with magnitudes and directions indicated by the color 
wheel in Figure 3 (bottom right). Information outside and inside parentheses characterizes currents in 
the equatorial (Figure 3, middle) and meridional plane (Figure 3, bottom), respectively. The maximum 
scale is 3 nA m -2 in the equatorial plane and 30 nA in the meridional plane. Note that diamag- 
netic currents are not represented in this diagnostic figure. 
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Figure 4. Simulated plasma densities in a cut through the Kwajalein latitudinal plane. Results are shown 
for simulation times of (top) 00:35, (middle) 01:10, and (bottom) 01:43, respectively, after a start time 
to = 09:30 UT (LT = UT + 1 1 h at the horizontal center of the simulation). 
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Figure 5. (top) Electron density (cm -3 ) and (bottom) zonal electric field (mV m 1 ) at 11:13 UT obtained 
from a longitudinal cut through the simulation at 379 km altitude over the Kwajalein sector (~400 km 
apex altitude). 


[ 25 ] The divergence in the vertical current density is sup- 
ported by a vertical electric field shown in Figure 3 (middle 
right). Above the bottomside, the field is directed down- 
ward, consistent with an eastward E x B drift in the same 
direction as the dynamo wind. Across the bottomside, the 
vertical electric field changes direction, producing a sheared 
zonal E x B drift, which provides free energy for the CSI 
below the shear node due to the differential velocity between 
plasma and neutrals. 

[ 26 ] Figure 3 (bottom) shows equatorward currents in both 
hemispheres on field lines that thread the bottomside in the 
altitude range from 250 km to 280 km. These currents feed 
the upward current driven by the F region zonal wind seen in 
Figure 3 (middle). Above the bottomside, the meridional 
current is distributed across all the field lines and is poleward 
in both hemispheres, serving to close the current loop driven 
by the F region dynamo. 

[ 27 ] Figure 4 shows the results for times matching, 
roughly, the three ALTAIR scans and C/NOFS passes. Note 
that the vertical axes indicate the apex altitude, i.e., the 
altitude mapped to the magnetic equator. For example, the 
magnetic field line that crosses over Kwajalein at 380 km 
altitude maps to 400 km altitude at the magnetic equator. 

[ 28 ] Figure 4 (top) depicts the evolution of ionospheric 
irregularities in the absence of a background electric field 
due entirely to the collisional shear instability. The spatial 
scales of the irregularities forming at the base of the bot- 
tomside range between 30 and 50 km. These irregularities 
cannot undergo vertical development and remain confined to 
altitudes where the plasma flow is retrograde. 

[ 29 ] By the time of Figure 4 (middle), the transient phase 
of the CSI was ending, and the asymptotic phase was 
underway. The 30-50 km waves at the base of the bottom- 
side have coalesced into large-scale waves with scale sizes 
between 100 and 200 km. Under the action of an increasing 
background zonal electric field, the bottomside irregularities 
begin to exhibit vertical development, more so on the eastern 


side of the simulation volume, which has experienced an 
eastward background electric field for more time. 

[ 30 ] In Figure 4 (bottom), under the influence of a signif- 
icant background electric field, the irregularities penetrate to 
the topside where gravity-driven currents can contribute 
further to growth. Some of the original intermediate-scale 
structure from Figure 4 (top) in the simulation survive, con- 
tributing to the overall fine structure. Bifurcation is evident, 
as well as secondary instabilities growing mainly on the 
western walls of the primaries. Some descending, wedge- 
shaped plasma enhancements accompany the ascending 
depletions. The characteristic westward tilt of the irregulari- 
ties is associated with the combination of vertical and zonal 
currents and the combined CSI and GRT instabilities. As the 
simulation was initialized with white noise, and as GRT is a 
universal instability (nearly scale size independent), the scale 
sizes of the irregularities evident here are a consequence 
of CSI. 

[ 31 ] Had a satellite passed through the simulated iono- 
sphere at 379 km altitude over the Kwajalein sector 
(''-400 km apex altitude) at about 11:13 UT, the electron 
density and electric field perturbations would have appeared 
as shown in Figure 5. Plasma depletions are present prefer- 
entially in the eastern half of the simulation grid. The west- 
ern walls of the irregularities have steeper density gradients 
than the eastern walls. Perturbations in the zonal electric 
fields of the order 3-5 mV m 1 are present, collocated with 
large and small plasma depletions alike. Since the pass is a 
few kilometers above the top of most of the plasma deple- 
tions, relative perturbations in the density are small in all but 
one case, but the fringing fields penetrate to higher altitudes 
and are detectable in all cases. 

6. Discussions 

[ 32 ] In simulation, the combination of CSI and GRT pro- 
duced growing waveforms with characteristics that matched 
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the observations made by ALTA1R radar. The dominant 
plasma irregularities with (initially) decakilometric-scale 
took approximately 100 min to evolve from the bottomside 
up to the F peak. By the time they reached the F peak, they 
have coalesced into waveforms with 100+ km scale sizes. 
This behavior was observed in ALTAIR scans and simula- 
tion results. 

[33] The transient response of CSI explained ~30-50 km 
wave observations, a feature commonly observed by back- 
scatter scatter radars in the form of patchy bottom-type ESF. 
These intermediate-scale irregularities were observed ini- 
tially in both the simulations and the radar observations but 
did not reappear after they coalesced into emergent 1 00+ km 
scale irregularities. This result agrees with observations 
made at other longitudinal sectors, e.g., at Jicamarca 
[Hvsell, 2000], 

[34] During the asymptotic growth phase, when 
intermediate-scale depletions dominated the scenario, some 
of the irregularities remained confined to the bottomside, 
while others experienced more vertical development. Had 
these irregularities been observed with coherent scatter, e.g., 
by the Jicamarca radar, they would have been labeled bot- 
tomside ESF. The irregularities in the simulation in Figure 4 
(middle) and radar data in Figure 1 (top) are comparable in 
terms of scale sizes and altitude spans. 

[35] The last stage of development of ESF occurred when 
the vertical extension of the plasma depletion was large 
enough that the irregularities crossed through the peak into 
the topside. Similar features are observed in both the simu- 
lation (Figure 4, bottom) and the radar data (Figure 1, bot- 
tom). Depletions that reached the F peak were spaced by 
hundreds of km and displayed tilts close to 30° westward 
from zenith. 

[36] Measurements made by IVM showed that electric 
field perturbations of the order of a few mV m -1 were 
sometimes observed even when depletions in the densities 
were not detected. Similar features were present in the sim- 
ulation and seemed to be a natural characteristic of the ESF 
morphology. The simultaneous measurements of small rela- 
tive perturbations in the density and large zonal electric field 
perturbations seems to be a consequence of a satellite pass 
near the top of a plasma depletion. Since the depletion itself 
was not intercepted, it was not observed, but the fringing 
fields penetrated to higher altitudes and were detected. 

[37] Zonal neutral winds in the bottomside and valley 
regions exerted control over initial irregularity formation. 
This sensitivity to the winds was expected, since the growth 
rate of the collisional shear instability depends on the 
counter streaming flow of the neutrals and plasma [Hysell 
et al., 2006]. As wind measurements in the lower thermo- 
sphere and mesosphere show considerable variability [Larsen, 
2002], sensitivity tests to winds at those altitudes need to be 
performed to evaluate their effects on ESF initiation. 

7. Summary 

[38] The combination of generalized Rayleigh-Taylor and 
collisional shear instability produced growing waveforms 
with characteristics matching satellite and radar observations 
in terms of their gross morphology and rates of develop- 
ment. Growing intermediate- (~30 km) and large-scale 
(100 + km) waves resulted from the nonlocal and nonnormal 


characteristics of the system and were not imposed a priori. 
The transient response of CSI explained rapid intermediate- 
scale wave emergence at the base of the F region, whereas 
the steady state behavior seemed to account for the large- 
scale waves that ultimately dominated. The three different 
types of ESF (bottom-type, bottomside, and topside) 
appeared to be a result of the time sequencing of the tran- 
sient and steady state response of CSI in relation to the zonal 
electric field. 

[39] A major objective of this research was to simulate 
the electrodynamics of ionospheric plasma irregularities, 
including those in ESF, in the direction parallel to the geo- 
magnetic field. The results of ESF simulation presented here 
suggests that a full 3-D treatment is necessary in order to 
track the distribution of parallel currents and the parallel 
closure of currents forced by dynamos. Parallel currents are 
found to have effects on the overall circulation and con- 
tribute to ESF morphology. Studies are underway to quan- 
tify the effects and mechanisms at work. 
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